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Abstract 

In this work, we investigate the existence of transition state resonances on atom-diatom reactive collisions from a time- 
dependent perspective, stressing the role of quantum trajectories as a tool to analyze this phenomenon. As it is shown, 
when one focusses on the quantum probability current density, new dynamical information about the reactive process 
can be extracted. In order to detect the effects of the different rotational populations and their dynamics/coherences, we 
have considered a reduced two-dimensional dynamics obtained from the evolution of a full three-dimensional quantum 
time-dependent wave packet associated with a particular angle. This reduction procedure provides us with information 
about the entanglement between the radial degrees of freedom (r, R) and the angular one (7), which can be considered 
as describing an environment. The combined approach here proposed has been applied to study the F+HD reaction, for 
which the FH+D product channel exhibits a resonance-mediated dynamics. 

Keywords: Transition state resonance, Bohmian mechanics, quantum trajectory, time domain, quantum 
hydrodynamics, survival probability 



1. Introduction 
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The analysis in the time domain of processes and phe- 
(— I nomena that appear in chemical reactions constitutes nowa- 
O j days a strong field of research both experimentally and 
l— ~~ 1 ^theoretically. From the experimental point of view, the 
femtosecond techniques developed in the 1990s have al- 
lowed us to monitor the process step by step, from re- 
00 actants to products, with an incredibly high resolution. 
This has lead to an impressive amount of theoretical work 
-based on wave packet (WP) propagation techniques. In 
. general, these studies are aimed at computing very long 
time series, from which energy/frequency spectra are ob- 
gained at a high resolution (the longer the timescale con- 
I sidered, the finer the features one can define in the corre- 
^ • sponding energy /frequency spectrum), which are later on 
, compared with the experiment. Within this scenario, little 
attention is eventually paid to the evolution of the WP it- 
?H self, apart from the frame-to-frame picture of the process. 



Not much information is usually extracted from the corre- 
sponding probability density, except for where the latter 
displays maxima or, on the contrary, presents regions of 
voids (nodal regions). 

In order to obtain an alternative description of the dy- 
namics, one can also analyze the quantum probability cur- 
rent density. Usually this quantity is only considered as 
a quantum flux to determine how much of the probability 
density flows into or out a certain region of the configu- 
ration space defining the system of interest. It also gives 
us an idea on how the former evolves in time, just as cur- 



rents or streams allow us to visualize how a river flows 
along its bed. This additional information has neverthe- 
less not been much exploited in the literature, except for 
the study of magnetic properties (because of the needs im- 
posed by the nature of the magnetic field, which is source- 
less and therefore has to be considered in hydrodynamic 
terms). In this sense, Bohmian mechanics P, a re- 
formulation of standard quantum mechanics in terms of 
trajectories or, equivalently, streamlines, results of much 
help. This approach enables us to take the analogy of 
the river further away: the trajectories behave like the 
paths pursued by some tracer particles that travel with 
the stream, thus giving us a precise description of the flow 
dynamics. This information cannot be obtained directly 
from the quantum probability current density. It provides 
us with an insight at a local level of the stream, but does 
not say anything about how a particular point-like particle 
evolves when it starts from a certain initial position in the 
corresponding configuration space. Within this approach, 
quantum trajectories (QTs) or streamlines are mainly ob- 
tained through any of two approaches Q: synthetic 
and analytic [6|,|jj- Within the former, QTs are computed 
after solving simultaneously the quantum Hamilton- Jacobi 
equation and the continuity one, while in the latter (used 
in this work) one starts from Schrodinger's equation and 
then the QTs are obtained from the phase of the wave 
function. 

The use of Bohmian mechanics with interpretational 
purposes can be traced back to 1926, when Madelung Q 
tackled the issue of the interpretation of the wave function 
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(in fashion at the time) by reformulating quantum me- 
chanics in hydrodynamic terms. These ideas underwent 
a rebirth in the 1970s through the works of Bialynicki- 
Birula d [lOl and Hirschfelder Within this hy- 

drodynamical framework, treating the probability density 
as a quantum fluid, the chemical reactivity of collinear re- 
actions was formerly studied by the end of the 1960s and 
beginning of the 1970s by McCullough and Wyatt [l5l-fl7|. 
This can be considered the starting point for the use of 
this theory in chemical physics, first at an analytic level 
and then at a synthetic one by Lopreore and Wyatt 0, [B[ . 
In this latter case, it constituted the source from which a 
series of numerical algorithms aimed at obtaining the evo- 
lution of quantum systems were developed, thus leading to 
the so-called QT methods 0] . In the same direction, more 
recently Sanz et al. have analyzed [HI, [l9[ the dynamics 
associated with a prototype of reactants-to-products re- 
action described by the so-called Miiller-Brown potential 
energy surface (PES) [20j. On the other hand, the QT 
methodology has also been exploited by different authors 
21 - 3l| to understand the magnetic properties of molecules 
within a framework that encompasses electronic structure 
and topology. 

Different theoretical methods have been developed in 
the literature to analyze, to explain and to understand the 
formation of transition state resonances. One of the rea- 
sons for such an interest relies on the crucial role played 
by these processes as intermediate stages in chemical reac- 
tion dynamics and aggregation processes to form atomic 
clusters. In this sense, the F+HD reaction constitutes a 
paradigmatic system. This system has been actively in- 
vestigated in the last decade due to the the observation 
through molecular beam experiments of a resonance for 



the reactive channel F+HD HF+D [3J,|33|. In the dif- 
ferential cross section, this dynamical feature manifests as 
a sharp variation of the peaks along the forward-backward 
scattering directions and, in the corresponding integral 
cross section, as a step-like profile. Quantum mechani- 
cal (QM) calculations on the Stark- Werner (SW) ab initio 
PES [34{ showed that those observations are consistent 
with a resonance-mediated dynamics of the reaction at a 
collision energy near 0.5 kcal-mol -1 (0.022 eV). In partic- 
ular, the existence of a resonance with (uj-h = 3, i>hd = 0) 
quantum numbers for the F-H and H-D stretching modes, 
respectively, located at the transition state region was the- 
oretically proved [32L |33j • On the contrary, the dynamics 
of the other product channel, F+HD — » FD+H, was found 
not to display evidences of similar features. 

The specific location of this resonant feature, just above 
the entrance channel, minimizes the number of partial 
waves contributing to the integral cross section. This also 
explains why the peak observed is capable to survive to 
the averaging introduced when all such partial waves are 
taken into account. In addition, the resonance lies just 
below the threshold for reaction predicted by classical ap- 
proaches, as it has been shown by means of quasi-classical 
trajectory calculations 35, l36j]. Although the maximum 



peak associated with the resonance was not reproduced in 
those studies, the rotational motion was found to play a 
significant role in the overall reaction mechanism. Thus, 
in order to provide a correct description of the dynam- 
ics of the F+HD — > FH+D reaction, a complete three- 
dimensional (3D) QM approach is required [32J. 

In this work we have carried out a fully converged QM 
study for a zero total angular momentum (J = 0) by means 
of a 3D time-dependent wave packet (TDWP) method. In 
an attempt to further investigate the interpretive capabili- 
ties of trajectory-based approaches for this kind of reactive 
processes, a QT method has been employed in combina- 
tion with the propagated WP. Thus, besides the 3D WP 
propagation, projections onto specific values of the angu- 
lar coordinate (7) have been considered to run an alterna- 
tive QT study of the process. The reason to consider the 
projection-based analysis presented here arises as a way to 
avoid the complexity of the full 3D QT dynamics at a first 
stage, attacking the study of the resonance by only con- 
sidering the QT propagation within the two-dimensional 
(2D) radial subspace. In spite of its limitations, this 2D dy- 
namical study is shown to be insightful when compared to 
the corresponding complete 3D TDWP calculations. More 
specifically, this Bohmian analysis allows us to determine 
how the energy /probability associated with a particular 
angle flows and can be related to the resonance. 

This work is organized as follows. In Section [2] the 
theoretical details of the TDWP and QT approaches are 
explained. Results are shown and discussed in Section [31 
Finally, the main conclusions extracted from this work are 
summarized in Section 2J 

2. Theoretical framework 

2.1. The WP procedure 

The TDWP method has been described before (37j . 
so here we will only describe the main details involved. 
Accordingly, consider the total Hamiltonian, H = T + V , 
is described in mass-scaled Jacobi coordinates. With this 
choice, the kinetic energy operator, T, for J = reads as 
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where 1 and j are the orbital angular and diatomic rota- 
tional angular momentum, respectively, and V(R, r, R • r) 
is the potential energy function, which in this work corre- 
sponds to the SW PES [34|. 

The WP time propagation can be formally expressed 

as 



|^(f))=e- i ( fi - fc ) , /^ i (0)). 



(2) 



Here, this propagation is performed using the symmetric 
split-operator formula |38| . 
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In Eq. ((2J, the initial-state selected WP, 4>i(0), is assumed 
to be a product state with the form 



&(0) = G(R)cl> vj (r)Y l j M (R,r) 



(4) 



where G(R) is the localized translational wave function, 
(/)„j(r) is the rovibrational wave function for the selected 
initial HD state, and 
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is a bipolar harmonic. The translational wave function 
G(R) is given by a Gaussian WP of width AR and average 
incident momentum Hko, i.e., 



G{R) 



_ e -(iJ-fl ) 2 /2Afl^-ifc fl 



(6) 



As it is explained in Ref. [37J, a discrete variable rep- 
resentation is employed to describe the radial R and r 
coordinates, while a basis set of bipolar harmonics is cho- 
sen for the angular degree of freedom, 7. For asymmetric 
reactions, as in the case considered here, 7 defines sepa- 
rate regions that correspond to two different product ar- 
rangements 37| . In order to avoid numerical reflections as 
the WP propagates along the different channels associated 
with the reaction, absorption of the corresponding outgo- 
ing probability fluxes is also considered. In this regard, 
for the absorbing potential —ie(R, r) = ien(R) — ie r (r) in 
Eq. ([2]) , we have used the expression proposed in Ref. [37| 
for the potential originally developed by Manolopoulos 
[39I ] . This is applied at r\ and R\, specific values at the 
end of the corresponding grids. 

The values for the parameters used in the TDWP cal- 
culation throughout this work, such as parameters for the 
initial WP (the center Rq of the Gaussian function G(R), 
the incident wave vector fen, and the width ARo), the time 
step (r), the total propagation time {ttot), the maximum 
value of diatom rotational states (jmax), the step for the 
radial grids (8r and SR) and the position of the absorbing 
potential for the radial coordinates (i?i and r\) are given 
in Table [U 

2.2. A brief account on Bohmian mechanics 

Consider a system of mass m in Cartesian coordinates. 
In Bohmian mechanics, the wave function 'J associated 
with this system provides us with dynamical information 
about it at any point on configuration space and any time. 
More specifically, this information is encoded in the phase 
of 'J, as infers from the transformation relation 
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where p and S are the probability density and phase of 
respectively, both being real-valued functions. From 
Eq. (J7]), Schrddinger's equation, 



(8) 



Table 1: Parameters used in the TDWP calculations carried out here 
(see text for details). 



Parameter 


Value 


R Q (bohr) 


9.6 


fco (bohr -1 ) 


0.42 


AR Q (bohr) 


0.09 


r(fs) 


0.016 


ttot (fs) 


165 


Jmax 


95 


SR (bohr) 


0.061 


8r (bohr) 


0.052 


i?i (bohr) 


11.0 


r\ (bohr) 


8.0 



is recast as a system of two real coupled equations, 
dp 
dt 

dS_ (VS) 2 
dt 

where 
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is the so-called quantum potential. Equation (|9]) is the 
continuity equation, which rules the ensemble dynamics of 
a swarm of trajectories with initial positions distributed 
according to po; Eq. (|10[) is the quantum Hamilton- Jacobi 
equation, which describes the phase field evolution ruling 
the motion of quantum particles through the motion equa- 
tion 

v = r=— . (12) 
m 

The coupling between Eqs. ^ and (fTUI) through Q (or, 
equivalently, p) is the reason why quantum (Bohmian) 
dynamics is very different from its classical counterpart, 
where both equations are only coupled through S. This 
therefore constitutes a very important difference between 
the typical calculations based on classical trajectories, com- 
monly used in chemical reactivity, although they cannot 
reproduce quantum features such as tunneling or interfer- 
ence, and simulations based on Bohmian trajectories. In 
other words, this coupling is precisely the way how the 
wave function guides the motion of the particle and hence 
allows QTs to display true quantum features. 

If instead of trajectories, we are more interested in 
looking at the quantum system dynamics as the evolution 
of a quantum fluid, i.e., in a quantum hydrodynamical 
view, the magnitudes of interest will be the probability 
density, p = , and the probability current density, 
J = pv = p(S7S/m), related through the continuity equa- 
tion (|9j), as 



(13) 



Hence, instead of talking about trajectories, one speaks 
about the conservation of the quantum flow, where the 
analog of Eq. (ITU1) for v is a quantum Euler or Navier- 
Stokes equation and the solutions of Eq. ([T2"j) are regarded 
as fluid streamlines rather than trajectories. These stream- 
lines are obtained by integrating 



v = r = 
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which is formally equivalent to Eq. (TT21 . Accordingly, 
these lines would follow the flow described by the quantum 
(probabilistic) fluid which describes the system. 

In this work, we are going to explore the 2D restricted 
dynamics associated with the subspace (R, r) describing 
the F+HD ->• FH+D/FD+H reaction. Assuming the last 
two terms of the kinetic operator given in Eq. (TTJ can be 
considered as a part of some effective potential operator, 
the motion equations describing the evolution of the cor- 
responding reduced QTs will be 
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where S 7 is the phase associated with the 2D time-depen- 
dent reduced wave function (j)j(R, r; t). This wave function 
is obtained by projecting of the total 3D wave function 
written in Eq. ([2]) onto a given angle 7. In principle, this 
projected wave function obeys a continuous, smooth evo- 
lution in the reduced subspace and, therefore, one also 
expects the corresponding reduced QTs to display a con- 
tinuous evolution. Note that this is somehow equivalent to 
describe the dynamics associated with an effective Hamil- 
tonian 
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Therefore, although this reduced dynamics does not pro- 
vide a full picture of the 3D system, as we shall see it 
results very insightful to understand the mechanism of the 
resonance process in a simplified manner. 

It is interesting to stress here that, due to its depen- 
dence on all the remaining values of the angular coordi- 
nate, the evolution under the Hamiltonian H 7 can be un- 
derstood as affected by a sort of environment coordinate. 
No reduction of the latter by tracing over the associated 
full 3D density matrix or by means of related phenomeno- 
logical reduced model ioL 41 1 is used. 



2.3. Quantum trajectories for a Gaussian WP 

One of the relevant aspects of the process concerns 
the preparation of the initial state, for it will determine 
the subsequent evolution of the chemical species. In or- 
der to render some preliminary light on this issue, here 
we are going to briefly analyze the evolution of our 2D 
projected WP. The wave function describing our system's 



state is given by the product of a Gaussian WP in the 
i?-coordinate and the ground state of the asymptotic po- 
tential along the r-coordinate. This means that, until the 
wave function is not well inside the interaction region, we 
will essentially have a wave function which is a free expand- 
ing and onwards propagating Gaussian with basically the 
same width along r. Thus, we can simplify the analysis by 
only considering the evolution of a Gaussian wave function 
evolving along the i?-coordinate for times lesser than the 
time ti n t at which the interaction is strong enough as to 
couple the two radial degrees of freedom and the (rota- 
tional) angular ones. In this sense, consider t <C ijnt- The 
evolution of ([5]), when normalized, is then approximately 
described by 



G{R,t) = A t e 
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where A t = {2naf)- 1 / 4 and AR = V2ct . This WP 
moves along the classical trajectory R c i — vot (for simplic- 
ity, we have assumed a zero value for the initial R c i,o = 
0) and its (complex) spreading with time goes as a t — 
<7o(l + iht/2m<TQ). The initial velocity by vq = hk^/fi = 
(p/fi) and the initial energy by E = fi, 2 fcg/2/i. The phase 
(Bohmian action) corresponding to (fl"8|) then reads as 
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being the time-dependent real spreading. 

The role of the quantum phase in the reaction dyna- 
mics is very important and becomes very apparent through 
the analysis of QTs. An analysis of the early stages of the 
reaction results very insightful in order to understand the 
further time-evolution. As said above, these stages are well 
described by the Gaussian WP (TT5)) . Thus, introducing 
(|19p into (|15p and then integrating in time, we find 



R q {t) = R c i(t) H R q ,o 



(21) 



where R q $ is the initial position of the corresponding tra- 
jectory. Notice that, in order to get this result, we have 
assumed that <f>y(R,r;t) can be considered separable in 
G(R, t) and another function depending on r at these early 
stages of the evolution, which also means that 



,(*) 



r <z,o, 



(22) 



i.e., at short times QTs only evolve along the .R-coordina- 
te, their position remaining essentially constant along the 
r-coordinate. 

In Eq. (f2"Tj) we clearly distinguish the two contributions 
ruling the time dependence of G(R). First, a classical drift 
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which makes any QT to move alongside the correspond- 
ing classical path. This could be therefore regarded as a 
Bohmian classicality criterion, to some extent related to 
Ehrenfest's theorem (see below). The second contribution 
is a quantum fluid drift, which comes from the own (quan- 
tum) nature of the evolution of 4> 7 and describes the (free) 
expansion of the WP with time. Accordingly, QTs will 
separate among themselves at a nonuniform accelerated 
rate or velocity, 



dR q 
~dt 
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(23) 



with r q = 2/icr§/fi.. That is, contrary to the classical de- 
scription, quantum-mechanically free expansion involves a 
stage of accelerated quantum motion. This means that the 
classical limit will imply keeping this term relatively small 
and, therefore, that QTs will be essentially parallel to the 
classical one, R c i, in agreement with Ehrenfest's theorem. 

By inspecting (|2"01) . two timescales become apparent 
(42| |. If t <C t 9 , the WP spreading is meaningless {a t ~ Co) 
and QTs are all essentially parallel to R c i (Ehrenfest's or 
classical criterion), since R q {t) « R q $ + R c i(t). In this 
case, the classical drift becomes the leading term. As t 
increases, if Eq. (|2~Tj) is expanded linearly in time, we start 
to observe the action of the quantum component of the 
velocity. This leads to an incipient accelerated motion, 
which at early stages can be described according to the 
familiar expression from classical mechanics 

1 



R q (t) « R q fl +v Q t+ - a q t 



(24) 



although the acceleration a q = R q ,o/ T q here depends on 
the initial position R q o — the larger the distance \R q fi — 
R c i,o\, the larger the effects due to this quantum accelera- 
tion. As time increases asymptotically (t 3> r q ), we find at 
approaches a linear dependence on time and the WP cur- 
vature becomes invariant under translations in time (this 
is the so-called Fraunhofer regime in optics, where phases 
depend linearly on coordinates). In this case, QTs go as 
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i.e., the asymptotic motion is again uniform, but with the 
important difference that the corresponding (constant) ve- 
locity having a component proportional to the initial po- 
sition of the trajectory. This extra QM component is very 
important: in the case of diffraction by periodic struc- 
tures, for example, it is precisely the one associated with 
the different diffraction channels [43[ or with the transfer 
of crystalline momentum [44| |45[ . 

In terms of the 2D problem, this means that at the 
first stages of its evolution, <fi~ will undergo a kind of hy- 
perbolic expansion along the incoordinate, independent 
of r (where the trajectories will remain essentially steady) 
and 7 (for any angle we will observe a similar behavior 
initially). Actually, in order to determine the leading dy- 
namics from the initial state (i.e., whether the spreading 



will dominate over the propagation or vice versa) , with the 
values given in Table Q] we find that the value of the ra- 
tio vq/v s — 2fco<7o .46!] is approximately 0.027. This means 
that the WP (fT5)) will propagate along R much slower than 
it spreads. Therefore, essentially at least half of it will 
move backwards and will not contribute to the resonance 
(in other words, the associated amount of rcactants will 
not be reactive). Once the onwards part of the WP (i.e., 
that spreading towards smaller values of R) starts reach- 
ing the transition state region, the approach considered 
here will not be valid any more (we will be in the regime 
t > ti n t). In this case, it is also possible that, due to re- 



flection or interference [18j, [19j , some additional portion of 



reactants will also become inactive, bouncing backwards. 

3. Numerical simulations 

3.1. The model 

In the (R, r) subspace, the SW PES i.34] employed here 
to study the F+HD reaction is strongly dependent on the 
Jacobi angular coordinate, 7, as can be seen in Fig. [TJ 
where a contour-plot of this PES at three different angles 
is displayed. 

In panel (a), for 7 = 0°, the entrance F+HD channel 
leads smoothly to the FD+H product channel which ex- 
hibits two separated potential minimum regions, a feature 
due to the possible atom permutations on the correspond- 
ing diatom. As 7 increases these two regions start to get 
closer and finally merge on an unique minimum well at 
7 = 90° (see panel (b)). If we keep increasing the an- 
gle, the PES starts to describe the other possible prod- 
uct channel, FH+D. Again, the presence of a significantly 
large barrier creates two unconnected regions, with the 
separation getting at its highest for 7 = 180° (see panel 
(c)). These projections of the PES at specific angles reveal 
that the connection of the separated regions observed for 
7 = 0° and 180° is in fact only possible via the 7 ~ 90° 
direction. That means that, whereas this pathway is cer- 
tainly possible for the 3D evolution of the WP, restrictions 
are expected for the reduced 2D dynamics of the QTs. 

As mentioned above, the reduced 2D projection is ex- 
tracted from the 3D WP at each time step St by fixing the 
angular coordinate, 7. It is worth noticing that the thus 
obtained evolution of the 2D WP onto the (R, r) subspace 
will not be equivalent at all (unless the angular degree 
of freedom is completely decoupled) to the evolution of a 
strictly 2D WP along any of the potential surfaces dis- 
played in Fig. [TJ On the contrary, due to the coupling to 
the angular coordinate, the 3D WP will lead the evolu- 
tion of its 2D projection, avoiding motions that one could 
observe otherwise. 

The initial conditions for the corresponding QTs will be 
chosen taking into account the initial 2D projected prob- 
ability density at a specific value 70 for the angle, 



P70 = I07o (0)r = \GiR)^j(r)Y^ M ( 7 o)| 



(26) 
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Figure 1: 2D PES (slices) corresponding to the three angles selected 
here: (a) 7 = 0°, (b) 7 = 90°, and (c) 7 = 180°. Colored contours 
cover a range of energies between the minimum energy (deep blue) , 
around —1.5 eV, and a maximum energy of 3 eV (red). 



Although this initial-condition density distribution can be 
approximated by a Gaussian function in both directions, 
we have implemented a Monte-Carlo-like distribution based 
on defining a grid on it and populating each cell accord- 
ing to the corresponding average weight. This procedures 
renders a fair distribution, as can be seen in Fig. [21 which 




9.7 9.8 
fi(bohr) 

Figure 2: Monte-Carlo-like sampling of the initial conditions for the 
QTs using the initial probability density. The red bigger dots mark 
the position of the tracer QTs (see text for details). 



consists of about 12,800 QTs when we choose a 100x100 
grid centered on p~ l0 and impose a maximum of 20 tra- 
jectories for the maximum value of p l0 . Apart from that, 
a less fine grid is also chosen to monitor the evolution of 
relatively large samples of QTs. The center of each cell of 
this grid is also chosen as the initial condition for a QT, 
namely the tracer QT, which will help us to understand 
the coarse grained evolution of the initial WP. In Fig. [21 
these tracer trajectories appear as red full circles. 

Once the initial conditions are set up, QTs are obtained 
by integrating Eqs. (|T5|) and (IT6l) with the value of the pro- 
jected 2D WP that we obtained from the propagation from 
t to t+r of the full 3D WP. It is clear that those trajectories 
approaching the limits of the entrance and exit channels 
may undergo the effects of the absorption. In principle, 
one expects that the phase of the guiding wave function 
will fade gradually and without changing much with re- 
spect to the unabsorbed wave function. In order to avoid 
eventual numerical problems, the evolution of those tra- 
jectories that approach the vicinity of the grid boundaries 
is stopped. That is, each time a QT reaches a maximum 
value for R and/or r, its evolution gets frozen. Note that 
this procedure is physically reasonable, since any trajec- 
tory moving backwards along the F+HD channel or for- 
ward along the FD+H or FH+D channels will already be 
moving asymptotically and, therefore, will not contribute 
anymore to the reaction dynamics. 

3.2. Analysis of WP dynamics and survival probabilities 

In order to monitor the formation of the transition 
state resonance and its survival, as well as the passage 
to products, when working within the time-dependent do- 
main the most appropriate tool is to define a survival or 
formation probability [47j |. 



Ps(«)= / P 10 {t)dRdr. 
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This quantity gives us the "amount" of probability inside 
a certain region £ as a function of time. Furthermore, 
its variation can be used as an indicator of the velocity 
of the reaction under study. In our case, we can consider 
the survival probability associated with reactants, prod- 
ucts and the transition state resonance. This means that 
three £ regions have to be defined. For example, looking 
at Fig. [IJa), the transition state resonance region, Sysij, 
can be defined as the region enclosing the "elbow" of the 
PES; the reactants region, as the region enclosing only 
the entrance F+HD channel; and the products region, Ep, 
as the region above £ts.r (which would include both exit 
channels). 

It is worth mentioning that this approach differs from 
the spectral quantization analysis performed by Skodje et 
al. [32| . In that work, the transition state was probed by 
explicitly setting the initial WP at this region. Here, the 
survival probabilities are defined in terms of a WP which 
initiates its evolution from the reactant arrangement. 

In Fig. [3] we have plotted the three different proba- 
bilities for the three values of 70 considered above (see 
dotted lines). As can be noticed, the two extreme cases 
(70 = 0° and 70 = 180°) present a certain resemblance: 
a fast decay of the reactants probability between 30 and 
50 fs and an equally fast increase of the products proba- 
bility, although this happens much earlier and faster along 
the FH+D channel (70 = 180°). The transition state res- 
onance probability, though, acquires relatively low values. 
In the case of the T-shape configuration (70 = 90°), the 
spatial limitations of the PES for the product channels (see 
Fig. Hlb)) make that all the probability transferred goes 
to the intermediate resonance region. 

The resemblance between the collinear configurations 
and their difference with respect to the T-shape one can 
also be seen if we compute the average kinetic energy as- 
sociated with 



"7o> 



h 2 ( d 2 d 2 \ 



!*»(«)>. (28) 



or the squared modulus of its time correlation function, 



c(t) = |<<Mo)|<M*))| 



(29) 



In Fig. 0] we have plotted these two quantities as a func- 
tion of time for the three angles considered. As can be 
seen, in panel (a) we find that K increases quickly and 
then decays more slowly and exponentially for 70 = 0° 
and 180°, while it decreases and becomes almost constant 
for 70 = 90°. Nonetheless, the maximum for 70 = 180° is 
less pronounced than for 70 = 0° , which is connected with 
the the different dynamics expected for the two possible 
product arrangements. The existence of a resonance for 
the F+HD FH+D channel leads the WP to stay longer 
at the spatially confined transition state region, while the 
unperturbed evolution associated with a direct mechanism 
that governs the dynamics for the F+HD — > FD+H pro- 
cess. This distinct behavior could account for the larger 




Figure 3: Reactants (black), transition-state (red) and products 
(green) probabilities as a function of time, obtained from the 2D 
reduced WP (dotted line), QTs (solid lines) and QTs with absorp- 
tion (dashed lines). From top to bottom, these probabilities refer to 
the three different angles considered here: (a) 70 =0°, (b) 70 = 90°, 
and (c) 70 = 180°. 



kinetic energy released for this latter channel. Additional 
support to this hypothesis will be discussed in the next 
sections. Regarding C(t), we find a similar trend for the 
two limiting cases, with an abrupt decay around 30 fs that 
does not appear for the T-shape configuration. 
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Figure 4: Average kinetic energy (a) and squared modulus of the 
time correlation function (b) associated with the reduced 2D wave 
function 7o as a function of time for the three the different angles 
considered here: 70 = 0° (black), 70 = 90° (red), and 70 = 180° 
(green). 



A similar analysis of survival probabilities can be now 
carried out in terms of QTs. In this case, the meaning of 
survival probability becomes closely related, as in classical 
statistical treatments, to the number of trajectories 
(each associated with a system configuration) that remain 
inside the corresponding region at a given time in relation 
to the total number of trajectories N considered in the 
statistics 14 711 . 



W(t) 



N ' 



(30) 



In principle, if the dimensionality of the trajectory-based 
calculation would be the same as the WP one, W would 
approach V(t) as N — > 00 (given an initial sampling ac- 
cording to po). Since this is not the case here for the two 
simulations performed, W will give us information about 
the 2D dynamics, which may differ importantly from the 
full 3D WP simulation. This is precisely what we ob- 
serve in Fig. [3] with solid line: a very small portion of the 
initial configurations sampled will contribute to the reso- 
nance and to products, since most of the probability will 
flow backwards (see next Section). Of course, as some 



probability is being absorbed, we can also compute W 
in terms of a variable total number of trajectories, which 
gives us the total amount of trajectories that are still in- 
side certain boundaries (see previous Section). However, 
the corresponding results (dashed lines in Fig. [3]) show that 
the probability mainly concentrates along the reactants re- 
gions (indeed, an important portion of it flows backwards, 
towards large R), with no significant contribution from 
the resonance-forming or the appearance of products pro- 
cesses. 

3.3. Analysis of quantum trajectories 

In the previous Section we have analyzed "ensemble" 
or statistical properties, without paying any attention to 
the way how the probability flows given a certain angular 
projection. Thus, let us consider now the evolution of the 
QTs, in particular the tracer QTs for each one of the angu- 
lar configurations given above. In this regard, it might be 
insightful first to consider a few snapshots of the evolution 
of the reduced 2D WPs for each angle. This is illustrated 
in Fig. [5] 

The evolution for the 70 = 0° case ((a) panels on the 
left of Fig. [5]) shows the WP extending on both branches 
with a rich nodal structure around the intermediate range. 
These features though do not seem to remain as the prop- 
agation time evolves, thus merely indicating the existence 
of the expected direct process for the F+HD — > FD+H 
case. For 70 = 90° ((b) middle panels of Fig. [5]), the re- 
quired bridge between the separated branches shown by 
the PES at some other angular directions, the WP is re- 
stricted to a very localized area at the intermediate region. 
After spending some time there, a significant component 
of the WP moves backward into the reactant channel. The 
situation at 70 = 180° ((c) panels on the right of Fig. [5]) 
clearly reveals the presence of a resonance at the tran- 
sition state in the FH+D product channel. As the WP 
reaches the intermediate region a stable 3-node structure 
is formed, governing the passage to both the upper and 
lower structures of the exit channels. 

In Fig.[B]we present the QT counterpart of Fig.[5J The 
first remarkable feature, as indicated above, is that more 
than a half of the initial conditions lead to backward mo- 
tions, this being the cause of the high probability in the 
reactants region. This justifies the fact of using filtering 
techniques in order to detect transition state resonances, 
for the corresponding probability is so small compared to 
the reactants probability that it would be difficult to de- 
tect them otherwise. Now, if we compare the three an- 
gular configurations, we readily notice that the number of 
tracer trajectories remaining in the resonance region (or 
in a neighborhood) is much larger in the case of 70 = 180° 
(see panel (c))than for 70 = 0° (panel (a)), where all the 
trajectories display a direct passage to products. This re- 
sult is consistent with the features of the resonance ob- 
served in the TDWP calculation for 70 = 180° discussed 
above (see Fig. [5]). Moreover, the behavior of the QTs in 
this particular region, where it is spatially confined, would 
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Figure 5: Snapshots illustrating the time-evolution of the reduced 2D WPs corresponding to the three different angles considered here: (a) 
70 = 0°, (b) 70 = 90°, and (c) 70 = 180°. From top to bottom, the times associated with each frame of the three cases considered are: 
(a) 13.2 fs, 39.5 fs, and 72.4 fs; (b) 6.6 fs, 26.3 fs, and 65.8 fs; (c) 6.6 fs, 39.5 fs, and 65.8 fs. 



also support a smaller amount of average kinetic energy as 
compared with the FD+H channel as commented in the 
previous section. 

On the other hand, in panel (b) we notice that, for the 
T-shape configuration (70 = 90°), trajectories concentrate 
for a time around two different regions, which essentially 
coincide with the separated areas of the FD+H and FH+D 
channels (see Fig. []}. Furthermore, it is also worth stress- 
ing the fact that, since the dynamics is restricted to a 2D 
plane, the QTs starting in the F+HD channel will only 
be able to just one of the available regions for the FD+H 
channel (see panel (a)) or the FH+D one (see panel (c)). 
The component of the fully 3D WP which explores the up- 
per branch of the PES contours does not manage to carry 
the corresponding 2D QTs. The trajectory seems to per- 



ceive it but not populate it. That is, it is a sort of "empty 
wave" In order to make active this part of the wave 
function, one should consider that initially there is some 
set of steady initial conditions resting in the associated 
channel. Then, as this projection of the 3D wave function 
would start emerging, those initial conditions would start 
to evolve accordingly. 

A complementary picture of the time-evolution of the 
QTs arises when we look at each degree of freedom sepa- 
rately as a function of time. The corresponding plots are 
presented in Figs. [5J and [HI respectively. Taking into 
account this viewpoint, first we notice that for 70 = 0° the 
i?-component behaves according to the prescription for a 
Gaussian WP at early stages of the evolution, as indicated 
in Section 12.31 while the r-component is basically steady 
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Figure 6: Tracer QTs corresponding to the three different angles 
considered here: (a) 70 = 0°, (b) 70 = 90°, and (c) 70 = 180°. 
Non-reactive trajectories are displayed with red color, while reactive 
ones in (a) and (c) are in blue, as well as those transiently passing 
to the second channel in (b); with black, tracer QTs that reach the 
transition state region and then are scattered back. In order to facili- 
tate the identification of the initial conditions leading to each type of 
motion, they have been plotted in the insets with the corresponding 
colors. 



(this behavior can also be observed for the other two an- 
gular configurations) . Then, we find the three well-defined 
sets of QTs: those moving backwards (red), those moving 
backwards after reaching the "elbow" of the reduced PES 



B 1 




Figure 7: Time-dependence of the R (a) and r (b) degrees of freedom 
of the tracer QTs corresponding to 70 = 0° . Colors are as in Fig. [6] 



(black), and those passing to products (blue). If we go to 
70 = 90° (see Fig. [5]), we readily notice the appearance of 
two groups of trajectories, one exploring the beginning of 
the upper region observed in both product channels (blue) 
and another one into the lower structure (black). How- 
ever, due to the constraint of 2D motion, none of them 
can indeed penetrate into any of the channels and the tra- 
jectories are bounced backwards. Finally, in the case of 
70 = 180° (see Fig. 0, we observe a behavior analogous 
to the one seen for 70 = 0°, although now only a smaller 
portion of the QTs passes to products, the remaining ones 
keeping trapped around the corresponding PES "elbow" 
(see Fig. [Ha)). 

4. Conclusions 

The dynamics of the F+HD reaction for J = has 
been investigated by means of a TDWP calculation com- 
bined with a QT approach, where the latter was based on 
computing the trajectories associated with 2D WPs arisen 
from the projections of the full 3D calculations at given 
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Figure 8: Time-dependence of the R (a) and r (b) degrees of freedom 
of the tracer QTs corresponding to 70 = 90°. Colors are as in Fig. [6] 



values of the angular coordinate. The evolution of the tra- 
jectories has been found to be limited when an explicit 
angular correlation is required. However, the distinct re- 
action mechanism for each product channel is successfully 
described. A clear indication of the resonance-mediated 
pathway for the FH forming arrangement thus becomes 
evident within the QT analysis, whereas the 2D evolution 
of the trajectories for the FD+H product channel is con- 
sistent with the observed direct reaction mechanism. 

As mentioned in the Introduction, in this work we have 
focused on the reduced radial subspace in order to avoid 
the complexity of the full 3D QT dynamics. For the same 
reason, we have chosen the case J = 0. This does not 
mean that the method is limited either theoretically or 
computationally. Rather it constitutes a firs step in the 
development of a more general QT methodology aimed 
at analyzing resonance formation processes. Notice that, 
by proceeding in this way, one gains an insight on the 
dynamics that can be later on used to understand in depth 
the full 3D dynamics. In this regard, if the 2D reduced QT 
dynamics provides us with interesting information about 
the flux of probability in these processes for a fixed angle, 
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Figure 9: Timc-dcpcndcncc of the R (a) and r (b) degrees of freedom 
of the tracer QTs corresponding to 70 = 180° . Colors are as in Fig. [6] 



the 3D dynamics will allow us to understand how this flux 
is transferred from among different rotational populations. 
Nevertheless, work in the direction of recreating the full 3D 
dynamics is currently under development. 
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